#########################
#                       #
#     Simulations       #
#                       #
#########################
rm()
memory.limit(18000)


library(plm)
library(knitr)
library(survey)
library(matrixStats)
library(nlWaldTest)

#Create table for results
table.FD<- data.frame(`FD_Beta1_Type1`="FD Beta1 Type1")

rownames(table.FD)[1] <- "T"

for (k in 5:100) {
  #First, define the parameters used by the data-generating process:
  rho <- 0.8
  sig1 <- 1 
  sig2 <- 1 
  beta1 = 0
  LT <- k*5
  T<-LT+50
  reps <- 100000
  
  set.seed(123)
  
  #The following function generates a synthetic dataset of desired dimension (T time points) and distribution parameters:
  
  generate <- function(T, rho, sig1, sig2) {
    y <- matrix(0, T)
    xI <- matrix(0, T)
    xS <- matrix(0, T)
    for (t in 1:T) {
      xxS <- if (t>1) xS[t-1] else 0
      xS[t] <- rho * xxS + rnorm(1, sd = sig1)
      xxI <- if (t>1) xI[t-1] else 0
      xI[t] <- xxI + rnorm(1, sd = sig2)
      y[t]<- xI[t] 
    }
    FDxS <- matrix(0, T)
    FDy <- matrix(0, T)
    FDxS[1] <- 0
    FDy[1] <- 0
    for (t in 2:T) {
      FDxS[t] <- xS[t] - xS[t-1]
      FDy[t] <- y[t] - y[t-1]
    }
    data.frame(t = seq(1,T-50),
               FDx1=c(FDxS[(51):(T)]),
               FDy = c(FDy[(51):(T)]))
  }
  
  #Now we generate datasets T time points 
  
  ds <- replicate(n = reps,
                  generate(T = T,
                           rho = rho,
                           sig1 = sig1,
                           sig2 = sig2),
                  simplify = FALSE)
  

  #Now we fit the FD model to the datasets and save the results
  set.seed(421)
  lms<-lapply(ds,
              function(d) {
                lm(FDy~FDx1,
                   data = d)
              })
  
  #Let's check the sampled lm parameters:
  true_param <- c( alpha0 = 0, beta1 = 0)
  est_param <- sapply(lms, coef)
  resid <- sweep(est_param, 1, true_param)
  est_confint <- sapply(lms, function(lm.model){confint(lm.model, parm=c(2))})
  power_beta1<-ifelse(est_confint[2,]<0 | est_confint[1,]>0 , 1, 0)
  power<-c(mean(mean(power_beta1)))
  
  
  table.FD[k-3,]<- data.frame(`FD_Beta1_Type1` = power[1])
  
  rownames(table.FD)[k-3] <- (k*5)
}
write.csv(table.FD, "FD--type1.csv")


